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Abstract 

The paper presents analytical or semi-analytical solutions for the formation 
and evolution of localized plastic zone in a uniaxially loaded bar with variable 
cross-sectional area. A variationally based formulation of explicit gradient 
plasticity with linear softening is used, and the ensuing jump conditions and 
boundary conditions are discussed. Three cases with different regularity of 
the stress distribution are considered, and the problem is converted to a 
dimensionless form. Relations linking the load level, size of the plastic zone, 
distribution of plastic strain and plastic elongation of the bar are derived and 
compared to another, previously analyzed gradient formulation. 

Keywords: plasticity, softening, localization, regularization, variational 
formulation 



1. One-Dimensional Softening Plasticity Model 

For many materials, the stress-strain diagrams characterizing their me- 
chanical behavior exhibit the so-called softening branches, with decreasing 
stress at increasing strain (and thus with a negative tangent stiffness). The 
physical origin of this intriguing phenomenon is in the initiation, propaga- 
tion and coalescence of defects such as microcracks or microvoids. Softening- 
induced localization of inelastic processes into narrow zones often acts as a 
precursor to failure. Proper modeling of the entire failure process requires 
an objective description of the localized process zone and its evolution. 

Perhaps the most popular class of inelastic material models is represented 
by the theory of (elasto-)plasticity. The present paper focuses on the lo- 
calization properties of softening plasticity models. To allow for analytical 



Preprint submitted to arxiv 



January 13, 2013 



solutions, all considerations are done in the one-dimensional context, re- 
ferring to the case of a straight bar under uniaxial loading as the typical 
paradigm. However, the analysis is nontrivial due to the fact that a variable 
cross-sectional area is considered, and a regularized formulation of softening 
plasticity is used. 

1.1. Classical Formulation 

In the small-strain range, classical elastoplasticity is based on the additive 
split of the total strain into the elastic part and the plastic part. The elastic 
strain is linked to the stress by Hooke's law, while the plastic strain can 
grow only if the stress level attains the yield limit, which is mathematically 
indicated by zero value of the yield function. The oriented direction of the 
plastic strain rate is specified by the flow rule and the evolution of the yield 
surface (set of plastic stress states in the stress space) is described by the 
hardening/softening law. For simplicity, we assume linear softening, i.e., 
linear dependence of the current yield stress on the cumulative plastic strain. 
Description of the stress-strain relation by a bilinear diagram is certainly a 
rough approximation, but it can reflect the main features of elastoplasticity 
with softening and serve as a prototype model, for which analytical solutions 
exist. 

In the one-dimensional setting, the basic equation can be summarized as 
follows: 

a = Eee = E{e-ep) (1) 

/(cr, k) = |o-|-cry(K) (2) 
Cry(K) = CTq + (3) 

Ep = A Sgno- (4) 
k = \ep\ (5) 
A>0, < 0, Xf{a,K) = (6) 

Here, a is the stress, e is the (total) strain, Se and Ep are its elastic and plastic 
parts, E is the elastic modulus, / is the yield function, ay is the current 
yield stress, do is the initial yield stress, H < is the softening modulus, A 
is the plastic multiplier and k, is the cumulative plastic strain. The overdot 
denotes differentiation with respec t to time. A more detailed d i scussi on of 



this specific problem is available in lJirasek. Zeman and Vondfeid (120101) and 



a broad backgroun d of th e theory of plasticity e.g. in iLublinerl (jl990l ) or 
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If we restrict attention to tensile loading (with possible elastic unloading, 
but never with a reversal of plastic flow), then the plastic strain e^, cumulative 
plastic strain k and plastic multiplier A are all equal. We will use k as the 
primary symbol for (cumulative) plastic strain and rewrite equations ([1]) and 
(ED as 

a = E{e- k) (7) 

K > 0, /(a, k) < 0, kfXa, k) = (8) 

The above equations refer to uniaxial tension, but formally the same 
framework can be used for the one-dimensional description of a shear prob- 
lem. Normal stress and strain are then replaced by shear stress and strain. 
Young's modulus E by the shear modulus G, and the tensile yield stress by 
the shear yield stress. 

1.2. Standard Gradient Formulation 

It is well known that softening is a destabilizing factor that may lead to lo- 
calization of dissipative processes (in our case of plastic yielding) into narrow 
zones. For classical continuum formulations with local dependence between 
stress and strain, the thickness of such localized process zones is undeter- 
mined and may become arbitrarily small. The undesired consequence is that 
the structural response becomes excessively brittle and numerical simulations 
suffer by pathological sensitivity to the discretization parameters such as the 
size of elements used by the finite element method. This has to be avoided, 
e.g. by introducing a regularization technique which enforces a nonzero thick- 
ness of the localized process zone and thus nonvanishing dissipation during 
the failure process. 

In the one-dimensional setting, negative plastic modulus H always leads 
to localization of plastic strain. Consider a straight bar with perfectly uni- 
form properties, subjected to uniaxial tension (induced by applied displace- 
ment at one bar end). The response remains uniform in the elastic range and 
also during plastic yielding with a positive plastic modulus. For a negative 
(or vanishing) plastic modulus, uniqueness of the solution is lost right at the 
onset of softening (or of perfectly plastic yielding). Stress distribution along 
the bar must still remain uniform due to the static equilibrium conditions (in 
the absence of body forces), but a given stress level can be attained by soft- 
ening with increasing plastic strain, or by elastic unloading with no plastic 
strain evolution. Which cross sections unload and which exhibit softening re- 
mains completely arbitrary, and there is no lower bound on the total length of 
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the softening region(s). Therefore, infinitely many solutions exist, including 
solutions with plastic strain evolution localized into extremely small regions. 
Even if the nonuniqueness of the solution is removed by a slight perturbation 
of the perfect uniformity of the bar, the problem with localization of soften- 
ing into arbitrarily small regions (in fact into the weakest cross section) still 
persists. 

Commonly used regularization techniques overcome the problem by suit- 
able enrichments of the governing equations. Such enrichments typically 
introduce at least one additional parameter with the dimension of length (or 
a parameter which can be combined with the traditional ones such that the 
result has the dimension of length). This parameter reflects the intrinsic 
length scale of the material and is related to the size and spacing of major 
heterogeneities in the microstructure. The size of the process zone is then 
controled by the choice of the length scale parameter. 

In principle it is possible to construct regularized models with enriched 
kinematic and equilibrium equations, e.g. strain-gradient plasticity or Cosserat- 
type models. From the practical point of view it is more convenient to limit 
the enrichments to the constitutive equations describing the material behav- 
ior and to keep the kinematic and equilibrium equations unchanged. This 
class of approaches is usually referred to as nonlocal continuum theories in 
the broad sense. Nonlocality of the stress-strain relation can be introduced by 
weighted spatial averaging of suitably chosen internal variables, or by incor- 
poration of gradients of such variables into the constitutive description. Here 
we focus on the latter case, in particular on its typical representative — the 
second-order e xplicit gradient model that evolved from the work of Aifantis 



and colleagues lAifantisI ( ll984j ) 



The explicit gradient formulation of elastoplasticity is based on incorpo- 
ration of a term proportional to the Laplacean of cumulative plastic strain 
into the softening law ([3]). In the one- dimensional setting, the Laplacean 
reduces to the second derivative and the enriched softening law reads 

ay(K) = (Jo + F(k + /\") (9) 

where / is a new parameter with the dimension of length. 

In a bar with perfectly uniform properties (cross section, yield stress, 
softening modulus, etc.) and in the absence of body forces and inertia forces, 
the stress is constant along the bar. The plastic zone, Ip, is characterized 
by growing plastic strain {k > 0) and vanishing value of the yield function 
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(/ = 0). Since the yield function is given by we conclude that the 
yield stress must be constant inside the plastic zone, and then leads to a 
second-order differential equation with constant coefficients and a constant 
right-hand side: 



Pk"{x) + k{x) 



H 



(10) 



As shown e.g. in de Borst and Miihlhaud ( 1992 ). the (most localized) plastic 
zone is an interval of length 27r/, arbitrarily placed along the bar. 

Analytical sol utions for several types of bars with variable cross sections 
were presented in iJirasek. Zeman and Vondfejd (120101 ) . The governing equa- 
tion 

F - aoA{x) 



I^A{x)k"{x) + A{x)k{x) 



(11) 



H 

was constructed from (|T0|) by setting cr{x) = F/A{x), where A is a function 
describing the distribution of the cross-sectional area along the bar, and F 
is the axial force transmitted by the bar, which is constant (independent of 
x) because of static equilibrium. 

In the present paper, we will use a modified formulation of the one- 
dimensional gradient plasticity model, constructed by a variational approach. 
Analytical or semi-analytical solutions will be derived and compared to the 
results for the standard gradient formulation based on ( fTTl) . An important 
advantage of the variational formulations is that it permits a consistent treat- 
ment of problems with discontinuous data, e.g. with a jump in the cross- 
sectional area (leading to a jump in the stress field). 

1.3. Variational Gradient Formulation 



itv model considered here is inspired bv the work of Miihlhaus and AifantisI 


(1991): Valanis (1996): 


Svedberg 


[(1996) 


: Svedbere and RunessonI (1997|.,1998il 


Polizzotto, Borino and Fuschi ( 


1998h: 


Borino, Fuschi and Polizzotto (119991): 


Liebe and Steinmann 


(2001). In the one-dimensional case, it is derived from 



the functional 



n(M,K) = [ lEA{u-Kydx+ [ \HA[k^ -fn^) dx + 
J c J c 



+ I AaoK dx — J Abu dx 



(12) 
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in which L denotes the interval that represents the entire bar, A is the cross 
sectional area and h is the prescribed body force density in the longitudinal 
direction (per unit volume), introduced just for the sake of generality but 
later set equal to zero. The first integral in ( !T2|) can be interpreted as the 
elastic strain energy, the second as the plastic part of free energy, the third 
as the dissipated energy and the fourth as the potential energy of external 
forces. 

Functional 11 is considered in the space of all sufficiently smooth displace- 
ment fields u that satisfy the geometric (essential) boundary conditions on 
the boundary 9£, and all sufficiently smooth and nonnegative plastic strain 
fields K. In formal mathematical language, the domain of definition of 
functional 11 is the space V = Vu x V,^ where 

Vu = {u & Hi{C) I u = u on dC in the sense of traces} (13) 
= G Hi{C) I K > almost everywhere} (14) 

This means that the functions describing the displacement and the plastic 
strain must be square-integrable and possess square-integrable generalized 
first derivatives, but continuity of the first derivatives and existence of the 
second derivatives are not apriori required. 

Due to the lack of convexity, the analysis can hardly rely on global mini- 
mization of functional 11. Nevertheless, it is reasonable to expect that stable 
solutions of the problem are associated with local minima of 11. The subse- 
quent derivations will be based on necessary conditions of a local minimum, 
in particular, on nonnegative values of the first variation (Frechet derivative) 
of functional 11 corresponding to all admissible variations of fields u and n. It 
will be demonstrated that such an approach leads to a consistent set of con- 
ditions that describe the problem and include the equilibrium equation, the 
complementarity conditions governing the plastic fiow, as well as appropriate 
boundary conditions at the physical boundary and regularity conditions at 
the elasto-plastic interfaces. A complete analysis should also pay attention 
to the second variation, which is related to stability issues. Analytical con- 
ditions for a non-negative second variation, derived for the simplest case of 
a bar with uniform properties, are presented in Appendix A. 

Strictly speaking, the variationa l approach sho uld be applied in an incre- 
mental fashion, as discussed e.g. by iPetryk fl2003[ ). However, for the present 



purpose it is fully sufficient to consider a total formulation. It turns out that, 
for one-dimensional problems with expanding or stationary plastic zones. 
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the parameterized solutions constructed in this way do not violate the irre- 
versibility constraints and thus represent physically admissible responses to 
given loading scenarios. 

The first variation of functional 11 defined in f|T2|) can be expressed as 



SIl{5u, 5k;u, k) = J EA{u' — k){5u' — 5k) dx + 

+ j HA{k5k-1^k'5k') dx + 

+ J Aao6K - y ^bSu dx (15) 

where 6u is the displacement variation (difference between two admissible 
displacement fields taken from Vu) and Sk is the variation of plastic strain 
(difference between two admissible plastic strain fields taken from V^). Inte- 
gration by parts of the terms with 6u' and 6k' leads to 

SU{Su,Sk;u,k) = - J [{EA{u - k))' + Ab]Su dx 

+ J [HAk + {HAl\y + Aao - EA{u' - k)] Sk dx 
+ EA{u' - K)n5u - ^[[EA(n' - k)5u]]^^ 

dC i 

- Y HAPkuSk + (16) 

dC i 

where Yldc stands for the boundary integral, in the one-dimensional setting 
reduced to the sum over two points at the boundary of the interval C, n is 
the "unit normal", equal to —1 at the left boundary and to 1 at the right 
boundary, and the sums with subscript i are taken over all points at which 
the quantity in the double square brackets has a discontinuity. The double 
brackets denote the jump of the quantity inside the brackets, defined as 

[[f]U= hm /(x)- hm fix) (17) 

X^XT X^X. 

As already explained, the first variation 511 must be nonnegative for all ad- 
missible variations bu and bK. By admissible variations we mean arbitrary 
changes of u and k for which n + 5n G K* and k + 5k G V^. 
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Variations 5u are arbitrary inside C but vanishing on the boundary. The 
expression multiplying 6u in the integral on the first line of f ll6p must be 
identically equal to zero, which provides the equilibrium conditions 

{EA{u' - k))' + Ah = Q (18) 

Of course, u' corresponds to the strain, u' — n is the elastic strain, E{u' — k) 
is the stress and EA{u' — k) is recognized as the axial force. Since 5u is 
arbitrary inside £, the jumps of EA{u' — k) must vanish, i.e., the axial force 
(not the stress) must remain continuous. 

The variation of k is not completely arbitrary, because n+5n must remain 
nonnegative. If we define the plastic zone Xp = {x G £ | n{x) > 0}, Sk 
can have an arbitrary sign inside Xp but must be nonnegative outside Xp. 
Therefore, the expression multiplying 6k in the integral on the second line of 
f lT6|) must vanish in Xp but outside Xp it is only constrained to be nonnegative. 
We recognize the resulting equation 

HAk + {HAPk')' + Aao = EA{u' - k) for x G (19) 

as the yield condition, and the resulting inequality 

RAk + {HAI^kJ + Aao > EA{u - k) ioi x e C\Xp (20) 

as the plastic admissibility condition. The advantage of the variational for- 
mulation is that the cases of variable sectional area A, softening modulus H 
or internal length / are covered in a systematic way, even in cases when some 
of these quantities exhibit discontinuities. From the last two lines of f|T6|) we 
obtain the corresponding jump conditions and also the boundary conditions. 

On the physical boundary dC, we get HAl'^n'n = if the boundary point 
belongs to Xp, or HAPn'ri < if this point does not belong to Xp. The 
first condition means that if the plastic zone contains a point of the physical 
boundary, the homogeneous Neumann condition k' = should be imposed at 
that point. The second condition means that if the boundary point remains 
elastic (k = at the boundary), the spatial derivative of plastic strain could, 
in principle, be nonzero. Since H < and Al"^ > 0, at the right end of the 
bar [n = 1) the derivative n' must not be negative. However, if k = at the 
right boundary and k > everywhere, k' cannot be strictly positive at the 
boundary and its only admissible value is again zero. 

The jump conditions imply that if the quantity HAI'^k'Sk is discontinuous 
at some point, its jump should be nonnegative for any admissible Sk. 
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• For points inside the plastic zone Xp, Sk can be positive as well as 
negative, and therefore HAt^n' must remain continuous. So in general 
we should not enforce continuous differentiability of plastic strain, but 
rather continuity of the product HAI'^k'. This is important when the 
spatial distribution of the softening modulus or of the sectional area is 
discontinuous. 

• For points outside the plastic zone, the variation 6k, is nonnegative 
and so the jump of HAI'^k' is in principle admissible but only if it is 
positive. Inside the elastic zone, k' vanishes and HAPk' has no jump 
at all. However, at the elasto-plastic interface (which is located at 
the boundary of the plastic zone), the derivative of plastic strain could 
exhibit a jump from zero value in the elastic zone to nonzero value in the 
plastic zone. At the left boundary of the plastic zone, the jump is equal 
to the value in the plastic zone. Since H is negative and At^ is positive, 
the limit of k' (as we approach the boundary of the plastic zone from 
inside) is allowed to be negative (or zero). However, k' < would mean 
that K < at some point inside the plastic zone (because k = at the 
boundary of that zone), which is not admissible. Similar arguments can 
be applied at the right boundary of the plastic zone, where the jump 
is minus the value in the plastic zone, and therefore k' is allowed to 
be positive (or zero), but a positive slope of the plastic strain profile is 
impossible to achieve without generating negative plastic strains inside 
the plastic zone near the boundary. So this discussion leads to the 
conclusion that the condition = should be imposed at the boundary 
of the plastic zone. 

In the absence of body forces, we set 6 = and equation (ITS]) implies that 
EA{u' — k) is constant along the bar (independent of the spatial coordinate 
x). Physically, this constant represents the axial force transmitted by the 
bar and therefore will be denoted as F. Equation f lT9|) is then rewritten as 

HA{x)k{x) + {HA{x)1'^k{x))' + A(x)ao = F for x G Xp (21) 

2. Bar With Piecewise Constant Stress Distribution 

Having presented the governing equations, we can proceed to localization 
analysis of a tensile bar with variable cross section. As the first case, consider 
a bar with piecewise constant sectional area (Fig. [T^). Suppose that the 
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bar contains a weak segment of length 2lg and sectional area A^, while the 
remaining parts have a larger sectional area A^/ (1 — /3) where /3 G [0, 1) is a 
dimensionless parameter. The origin of the spatial coordinate system can be 
placed into the center of the weak segment. The solution is then expected to 
exhibit symmetry with respect to the origin. 

Let us emphasize that the present analysis is strictly focused on one- 
dimensional modeling. Therefore, the stress distribution across each section 
is considered as uniform. Of course, for a real three-dimensional body con- 
taining notches, the stress field would have a singularity at the notch tip and 
the stress distribution across sections near that singularity would be highly 
nonuniform. However, we use the case of variable cross section as a paradig- 
matic example of the localization properties of a gradient plasticity model in 
cases with non-smooth and sometimes even discontinuous data. Therefore, 
the stress is expressed simply as the normal force divided by the sectional 
area. 

For the bar with a weak segment of length 21 g, the stress distribution is 
described by 



and has discontinuities at sections x = ±lg; see Fig. [2l As long as the stress 
remains below the yield limit, the response is purely elastic. The onset of 
yielding can be expected when the yield limit is attained, which happens at 
the elastic limit force Fq = A^cfq. For a softening plasticity model without 
any regularization, plastic yielding would localize into one cross section, the 
dissipation would vanish (because the plastic zone has zero volume) and 
the response would be extremely brittle. Regularization by the additional 
gradient term leads to a finite size of the plastic zone. 

2.1. Plastic Zone Contained in Weak Segment 

Let us first assume that the plastic zone Xp is fully contained in the weak 
segment, i.e., Xp C {—Igjg)- The yield condition ( 12T]) with A = A^ and 
F = Af-Oc can be rewritten as 




for |a;| < Ig 
for |a;| > Ig 



(22) 




(23) 
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(b) 



Figure 1: Tensile bars with (a) discontinuous distribution of stress, (b) continuous distri- 
bution of stress, (c) smooth distribution of stress 



of 



Oc (l-p)CTc 



la Ic 



"1 r 



J L 



Figure 2: Bar with a weak segment and the corresponding stress distribution 



This is a linear second-order differential equation with constant coefficients 
and a constant right-hand side, and its general solution is 

n{x) = + Ci cos y + C2 sin ^ (24) 

where Ci and C2 are arbitrary constants. Let Lp denote the length of the 
plastic zone and suppose that the plastic zone is centered at the origin, 
i.e., Zp = {—Lp/2, Lp/2). If this was not the case, the origin could simply 
be shifted to the center of the plastic zone. As explained at the end of 
Section 11.31 the plastic strain k must remain continuous and its derivative 
must vanish at the elasto-plastic interface, dXp, which consists of two points, 
X = ±Lp/2. Conditions K{-Lp/2) = 0, K{Lp/2) = 0, K,'{-Lp/2) = and 
K,'{Lp/2) = lead to C2 = 0,Ci = [a^-a^) /{H cos{Lp/2l)) and sin(Lp/2/) = 
0. The last condition means that the length of the plastic zone, Lp, must 
be an integer multiple of 2ttI. The shortest positive length of plastic zone 
is obtained if Lp = 2111. However, such a solution is admissible only if the 
plastic zone is indeed fully contained in the weak segment of length 2lg, i.e., 
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if Ig > Til. In that case, the solution is given by 

k{x) = I + t) ^ ^ = (25) 

[ foT X ^Xp 

This is the classical solution that describes localizat ion in a bar with perfectly 
uniform properties de Borst and Miihlhaus ( 1992I ). Analysis of the second 



variation of functional 11, presented in detail in Appendix A, reveals that 
this solution corresponds to a local minimum of 11. 

If the weak segment is sufficiently long with respect to the characteristic 
length of the material (longer than 2-7r/), the plastic zone will form inside 
that segment and the solution will not be affected by stronger parts of the 
bar. However, if the weak segment is shorter than the plastic zone in a 
perfectly uniform bar, solution (^5^ is not admissible and the derivation must 
be modified. 

It is convenient and elegant to convert the problem to a dimensionless 
format. For this purpose, we introduce the dimensionless spatial coordinate, 
C, = x/l, normalized plastic strain, Kn = —Hn/aQ, and dimensionless stress 
= (Jc/cTo. Note that 0, defined as the ratio between the stress in the weak 
segment and the yield stress, is at the same time the ratio between the axial 
force, F, and its elastic limit value, Fq, and thus will be referred to as the 
load parameter. The distribution of plastic strain in a uniform bar, given by 
f l25p . can be described in terms of the dimensionless quantities as 

UO = l 1'"'^^^' + '°'^^ (26) 
\ for ^ ^ (-7r,7r) ^ ^ 

This solution is valid for a bar with a weak segment of length 2lg provided 
that 2lg > 27r/, i.e., > vr where \g = Ig/l is a dimensionless parameter 
describing the ratio between the "geometric" characteristic length, Z^, and 
the material characteristic length, /. 

2.2. Plastic Zone Extending to Strong Segments 

Let us proceed to the case when 2lg < 27il, i.e., \g < re. Equation ( 123|1 
is now valid only for \x\ < Ig. In terms of the dimensionless quantities, we 
rewrite it as 

<(0 + ^^n{^) = 1-0, lel < A, (27) 

For simplicity, the derivatives with respect to ^ will still be denoted by primes, 
despite possible confusion with derivatives with respect to x. For the parts 
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of the plastic zone surrounding the weak segment, a similar equation with a 
modified right-hand side can be derived from 

<(O + '^n(O = l-0 + /30, Ag<|e|<Ap (28) 

Here, \p = Lp/2l is a dimensionless parameter characterizing the ratio be- 
tween one half of the plastic zone length, Lp/2, and the material characteristic 
length, /. Since the solution is again expected to be symmetric with respect 
to the origin, we construct it only in the weak segment and in one of the 
stronger segments adjacent to it: 

f 1 -0 + Cicos^ + Casing for - < ^ < 

f^niO = < (29) 
[ 1 - + + C3 cos{ + C4 sin^ for A^ < ^ < Ap 

Integration constant C2 must vanish due to symmetry. Constants Ci, C3 and 
C4 and the dimensionless plastic zone size Xp can be determined from four 
conditions: continuity of k„ and of Ak'^ at ^ = A^ and at ^ = Ap. These 
conditions lead to the set of four equations 

Ci cos \g — C3 cos Xg — C4 sin \g = f3(j) (30) 

-Ci(l -/3)sinAg + C3sinAg -C4COSA3 = (31) 

C3 cos Ap + C4 sin Ap = (1-/3)0-1 (32) 

— C3 sin Ap + C4 cos Ap = (33) 

which are linear in terms of Ci, C3 and C4 but nonlinear in terms of Ap. 
Since the load parameter, 0, enters the equations in a linear fashion, it is of 
advantage to reformulate the problem and solve for the integration constants 
and (f) in terms of Ap. Another parameter that affects the solution is A^, and 
so we can parameterize the solution by A^ and Ap and write 

C.iXg,Xp) ^ (34) 
C3(A„Ap) = i^-m"^^^^ (35) 
C4(A„Ap) = (36) 

'^^^^'^^^ = /sin A, 

^+D(A„Ap) 
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where 



D{\g,\p) = (1 — /3 sin^ Ag) sin Ap — /3 sin Ag COS COS Ap = 

= (l-|/3)sinAp + i/3sin(Ap-2A3) (38) 

Recall that parameter 0, defined as the ratio ctc/ coj can also be interpreted 
as the ratio between the axial force transmitted by the bar, F = Acac, 
and its limit elastic value, Fq = AcCTq. For a fixed A^, the dimensionless 
size of the plastic zone Ap can be varied as a parameter describing various 
stages of plastic zone evolution. The range in which the solution makes 
physical sense can be determined from the conditions that the initial value 
of the load parameter at the onset of yielding is = 1 and that = 
corresponds to complete failure. The relation between the load parameter 
and the dimensionless plastic zone size Ap is graphically illustrated in Fig. [3] 
for /3 = 0.5 and for several values of parameter Xg ranging from 0.01 to 2.5. 
The interesting range of Xg is between zero and tt because for Ag > vr the 
weak zone is long enough to allow formation of the complete plastic zone in 
its interior, and the solution is the same as for a bar with a perfectly uniform 
section Ac. 

When the axial force reaches the limit elastic value (i.e., when = 1), 
the entire weak segment of length 2lg starts yielding and the initial value 
of the dimensionless plastic zone size is Ap = Xg. Subsequently, the plastic 
zone expands continuously, and this initially happens at an increasing load 
level, provided that A^ < 7r/2, i.e., that the length of the weak segment 2lg 
is smaller than irl. The maximum force 

Fma. = ^ ^sinA, ^^^^ 



1 - (3(2 - (3) sin' A^ 
is attained when the dimensionless size of the plastic zone is 

1 — (3 sin^ A„ 

Ap,peafc = TT - arctau . — (40) 

P Sm Xg cos Xg 

After that, the axial force decreases to zero as the size of the plastic zone 
approaches the limit 

If the size of the weak segment exceeds vr/, the global response is softening 
right from the onset of plastic yielding. 
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' ' ' ' ' ' ^ — ^ 

0.5 1 1.5 2 2.5 3 3.5 

dimensionless plastic zone size, 

Figure 3: Relation between load parameter and plastic zone size for a piecewise constant 
stress distribution (bar with a weak segment) 




dimensionless coordinate, x/l dimensionless coordinate, x/l 

Figure 4: (a) Early stage, (b) late stage of evolution of plastic strain profile for a piecewise 
constant stress distribution 
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Substituting from (151|) - (157l) into we obtain the distribution of plastic 
strain parameterized by Ap, and integrating over the plastic zone we get the 
plastic elongation. An example of the evolution of plastic strain profile, 
calculated for /3 = 0.5 and Xg = 0.5, is shown in Fig. HI The plastic zone 
expands and, at each section, the plastic strain grows monotonically. 

It is also of interest to construct the load-displacement diagram for the 
entire bar. The elastic elongation is proportional to the axial force and the 
proportionality factor (bar compliance) depends on the bar length. There- 
fore, it is convenient to consider the contribution of plastic strain to the 
elongation separately, since the bar length does not affect it (provided that 
the bar is sufficiently long such that the full plastic zone can develop). The 
plastic part of bar elongation, 

/Lp/2 
k{x) dx (42) 
■Lp/2 

can be computed by integrating the plastic strain along the plastic zone. 
In the context of dimensionless description, it is natural to deal with the 
dimensionless plastic elongation 

/ Kn{x) dC = / = f = (43) 

where Kf = —(Tq/H is a material parameter that corresponds to the plastic 
strain at complete failure if the gradient terms are ignored. 

The plastic part of the load-displacement diagram, obtained by plotting 
the dimensionless load parameter against the dimensionless plastic elonga- 
tion Up/ Kfl, is shown in Fig. for fixed (3 = 0.5 and different dimensionless 
sizes of the weak segment, A^, and in Fig. \5jp for fixed Xg = 0.5 and different 
values of (3. These graphs explain the infiuence of both parameters on the 
shape of the load-displacement diagram. The weak segment can be consid- 
ered as an imperfection in a uniform bar, parameter A^ refers to the length 
of that imperfection and f3 to its "magnitude" (in the sense that larger /3 
means a more dramatic reduction of the sectional area). For short imperfec- 
tions, the initial structural hardening is very steep and the maximum axial 
force is close to the value calculated for a perfect bar. For somewhat longer 
imperfections, the structural hardening is more gradual and the maximum 
force is between the values that would correspond to the weaker and stronger 
sections, and for imperfection lengths exceeding nl (which is exactly one half 
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normalized plastic elongation, Up/lK, normalized plastic elongation, Up/lK, 



Figure 5: Plastic part of load-displacement diagram for a piecewise constant stress distri- 
bution for (a) different values of relative imperfection size Xg = Ig/l, (b) different values 
of imperfection "magnitude" /3 

of the plastic zone that would form in a perfectly uniform bar) the maximum 
load is dictated exclusively by the weak section and the load-displacement 
diagram is softening from the onset of yielding. 

3. Bar With Piecewise Linear Stress Distribution 

3.1. General Solution 

Now consider a bar with variable sectional area described by a function 
that is continuous but not continuously differentiable at the weakest section; 
see Fig. [TJ). To allow for analytical solutions, we define the specific distribu- 
tion of sectional area such that the corresponding stress distribution becomes 
piecewise linear. This is achieved by setting 

A{x) = (44) 

where lg is a parameter that sets the geometric length scale of the problem. 
Substituting (l44l) into (12T1) and dividing by HA{x), we obtain 

I | /t a;) + KX = "^^^^ ' '^ --^ fora;GXp-{0} 45 

lg — \X\ tilg n 

Equation psj) must be satisfied at all points inside the plastic zone, with the 
exception of point x = at which the sectional area is not differentiable. At 
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that point, conditions of continuity of k and An' have to be imposed. Since 
A is continuous, the latter condition actually means continuity of k' . After 
conversion to dimensionless form, the governing equation reads 

This is a second-order differential equation, but in contrast to ( 127|) or (128|) . it 
contains on the left-hand side an additional term with the first-order deriva- 
tive which has a non-constant coefficient, and the right-hand side is not 
constant. Still, an analytical solution in terms of special functions can be 
constructed. The derivation is presented in detail in Appendix B. The re- 
sulting general solution has the form 

= (A»-|fl)CfJi(A,-|?|) + (A,-|?|)C,±Y,(A,-|?|) + 

+1-^^^%^H,(A,-|«|) fo^i^I, (47) 

where Cf and are integration constants, Ji and Yi are the Bessel func- 
tions of the fi rst and second kind, resp., and Hi is the Struve function 



KorenevI fl2002t ) 



3.2. Particular Solution 

For simplicity, we have written (l47j) in a compact form, but the integration 
constants have different values in the "positive part" of the plastic zone, 
where ^ > 0, and in the "negative part" of the plastic zone, where ^ < 0. 
Therefore, we deal with four integration constants, C^, C^^, and C^, 
and with two additional unknowns that determine the position of the right 
and left boundary of the plastic zone. This makes a total of six unknowns, 
which can be determined from appropriate jump and boundary conditions. 
According to the foregoing analysis, the solution must remain continuously 
differentiable at the origin and at the two boundary points, which leads to 
six equations for the six unknowns. 

Due to symmetry, the problem can be simplified and it is sufficient to 
restrict attention to the positive part of the plastic zone, X+ = (0, Ap). The 
three relevant unknowns, , C2 and Ap, can be determined from three 
conditions, 

<,(0) = (48) 
«:n(Ap) = (49) 
<(Ap) = (50) 
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Substituting the general solution (H7|) into we obtain a set of three 

equations, 



^[-Yo(A, 


-A,)Ho(A,)+Yo(A,)Ho(A,-A,)] 


(54) 






7r[Jo(A9 - 


Ap)Ho(Ag) - Jo(A3)Ho(A3 - Ap)] 


(55) 




D+(A„Ap) 


2\g [Jo(Ag 


- Ap)Yo(Ag) - Jo(Ag)Yo(Ag - Ap)] 


(56) 




D+(A„Ap) 



-2AjC+Jo(A,) + C2+Yo(A3)] +07rHo(Ag) = (51) 
C+Ji(A,-Ap) + C+Yi(A,-Ap)-^Hi(A,-Ap) + -^ = (52) 

Z,Ag Ag — Ap 

-2A, [C+Jo(A,-Ap) + C2+Yo(A3-Ap)] +07rHo(A,-Ap) = (53) 

which are linear in terms of and but nonlinear in terms of Ap. Again, 
it is of advantage to reformulate the problem and solve for C^, and in 
terms of Ap and A^: 

Cl'(Ag,Ap) = 

C2^(Ag,Ap) = 

0(Ag,Ap) = 

where 

D+(A„Ap) = 2Ho(A,) + 7r(A,-Ap){[-Ji(A3-Ap)Yo(A3)+ (57) 
+Jo(A,)Yi(Ag-Ap)]Ho(A,-Ap) + 
+ [Jo(A, - Ap)Yo(A,) - Jo(A,)Yo(A, - Ap)] Hi(A, - Ap)} 

From symmetry with respect to the origin, we obtain integration constants 
corresponding to the negative part of the plastic zone as Cf = and 
C2 = C2 ■ The assumption of symmetry may seem to be restrictive, but 
it can be verified by (tedious) analysis of the complete problem with six 
equations and six unknowns that no nonsymmetric admissible solution exists. 

3.3. Results and Discussion 

Based on the solution constructed in the previous subsection, the evo- 
lution of the plastic zone can be analyzed, and the corresponding load- 
displacement diagram can be constructed. The plastic part of the load- 
displacement diagram is shown in Fig. |6l where the dimensionless load pa- 
rameter is plotted against the dimensionless plastic elongation, Up/lnf. 
Recall that (p = Of-ja^ = F/Fq is the ratio between the axial force F and 
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its value at the onset of yielding, Fq = A^aQ. The plastic elongation, Up, is 
normalized by a reference value luf, which would correspond to the plastic 
elongation of a bar of length I at complete failure if the solution remained 
uniform. 




normalized plastic elongation 



Figure 6: Plastic part of normalized load-displacement diagram for piecewise linear stress 
distribution 



During the elastic stage of loading, the axial force increases from to Fq 
and no plastic strain evolves. Thus the initial part of the diagram in Fig. |6l 
up to = 1, is vertical. For a bar with a uniform section, the continuation 
of that diagram would be a straight line that corresponds to linear softening, 
and the plastic elongation at complete failure (0 = 0) would be 2'KlKf. This 
is in fact a limit case of the present solution with \g oo. For finite Ag, 
i.e., for a bar with a variable section, the load parameter first increases and 
only later decreases. Complete failure is attained at larger elongations than 
for the uniform bar. This is represented in Fig. |6] by the solid curves, which 
correspond to different values of parameter = 3.2, 5, 10 and 50 (from 
top to bottom). Lower values of A^ correspond to stronger variation of the 
secional area and lead to higher peak loads an d higher elongations at failure. 
For co mparison, the solutions constructed in iJirasek. Zeman and Vondfejc 
( 120101 ) using a "standard" formulation, with the second term in ( 12T|) replaced 
by HAI^k", are plotted by the dotted curves. The present, variationally 
based formulation leads to qualitatively similar solutions, but with higher 
peak loads and higher elongations at complete failure. Note that, for the 
standard formulation, the elongation at failure is always the same as for a 
uniform bar, independently of the value of parameter A^. 
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Figure 7: Evolution of plastic zone size for piecewise linear stress distribution 



The evolution of the plastic zone size is illustrated in Fig. [7] by plotting 
the load parameter (f> against the dimensionless plastic zone size Xp — 
again for Xg = 3.2, 5, 10 and 50. For a uniform bar (Ac, — oo), the plastic zone 
would attain its full size Lp = 2ttI (corresponding to Xp = vr) immediately at 
the onset of yielding. In contrast to that, for the bar with a variable section, 
the plastic zone evolves continuously from the weakest section up to its full 
size, attained at complete failure. The fact that the plastic zone expands 
monotonically is important, because it justifies our tacit assumption that 
the material outside the current plastic zone has not experienced any plastic 
straining before. The full size of the plastic zone turns out to be somewhat 
smaller than fo r a uniform bar, which is diffe r ent fr om the standard solution 
constructed in iJirasek. Zeman and Vondfejc fl2010f ). plotted for comparison 
by dotted curves. 

The evolution of the plastic strain profile is plotted in Fig. |H]for Xg = 5. 
Again, solid curves represent the variational formulation and dot ted curves 
the standard formulation from lJirasek. Zeman and Vondrejc fl2010h . At early 
stages of plastic strain evolution, both formulations give almost the same 
result, but later the differences grow. 



4. Bar With Smooth Stress Distribution 

As the most regular case, we consider a smooth distribution of sectional 
area, with continuous derivatives of an arbitrary order. The origin of the spa- 
tial coordinate is again placed at the weakest section, where plastic yielding 
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dimensionless coordinate 



Figure 8: Evolution of plastic strain profile for piecewise linear stress distribution 



is expected to start. In Ijirasek. Zeman and Vondfejd ( l2010h . the function 
describing the sectional area was selected such that the resulting stress dis- 
tribution was quadratic. For the present purpose, it turns out to be more 
convenient to consider another special case in which the stress distribution 
is given by a Gaussian function, 



CTr.e 



(58) 



The sectional area is thus 



A{x) = A^e' 



'in 



(59) 



Since, for this case, the solution based on the standard formulation has not 
been published yet, let us present it before turning attention to the variational 
formulation. 



Ji-.l. Solution Based on Standard Formulation 



According to the standard formulation used in lJirasek. Zeman and Vondfejc 



( l2010l ). the yield condition is considered in the form 

HAfK"{x) + HAk{x) + Aao = F for xeXp (60) 

which differs from ( 12T1) by the second term. The corresponding differential 
equation for the plastic strain can be converted into the dimensionless form 

<(O + «:n(O = l-0e-«'/2^« (61) 
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The general solution of the corresponding homogeneous equation is a linear 
combination of cos^ and sin^, and the particular solution for the given right- 
hand side can be constructed by variation of constants. For the particular 
solution represented by 



UO=Ci{0 cose + C2(e) sine 
we obtain a set of two equations 



C-UO cose + 6-^(0 sine = 

-cue) sine + c^(e) cose = i 



-e/2\i 



from which 



c'i(e) 



cos 



e + ^y e-«'/2A?sin^ de 



sine — 



-e/2xi 



cose de 



(62) 



(63) 
(64) 



(65) 
(66) 



The integrals can be conveniently expressed by switching to complex func- 
tions. As shown in Appendix C, the resulting general solution of equation 
flM]) can be expressed as 



'^n(e) = Cicose+C2sine+l-0Age-« ^ 



(67) 

where Ci and C2 are arbitrary constants, i denotes the imaginary unit, and 



Fix) 



rat 



(68) 



is the so-called Dawson function; see e.g. lOlverl ( 1l997l ). 

In general, one should impose two boundary conditions (of vanishing value 
and vanishing derivative) at each boundary point and solve for two integra- 
tion constants Ci and C2 and two unknown coordinates of the boundary 
points. Due to symmetry, the problem can be simplified. Integration con- 
stant C2 must vanish, and the plastic zone Xp = (— Ap, Xp) is characterized by 
one unknown parameter, Ap. Boundary conditions K„(Ap) = and n'n{Xp) = 
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lead to a set of two equations, 



A. 



Ci cosA„ + 0-^e-^?/2A: 
V2 



A3 - «Ap 



v^A 



?A --X 



Ci sin Ap + 0^= e 
which can be written as 



g 



+ F 
- F 



CpCi + EFr^ 
SpCi + EFj^ 



A^ + iXp 
KV2 



1 




where c„ 



cos Ap, Sp 



sm 

1 
2 



Ap, E = v^A^e-^p/^A^ g^^^ 



A; + ^Ap 
KV2 



Ag + iAp 



+ F 



Ag - ^A, 
X,V2 



X,V2 



iXr 



X,V2 



(69) 
(70) 



(71) 
(72) 



(73) 
(74) 



are the real and imaginary parts of F (^"^^v^j- Equation ( ITOl) follows from 

the condition K'„(Ap) = with the derivative of F expressed as F'(a;) = 
1 — 2xF(x), which easily follows from the definition of Dawson function 
Solving for Ci and </>, we obtain 



Ci(Ap, Xg) 

0(Ap, Xg) 



' AgV2 



XnV2 



FrSp + FjCp 



V2 



iApp 



A„V2 



9 

i sin ApC'^'p/^'^a 



(75) 



E(F/jSp + FjCp) Xg g^App ^l+^p 



— e 



(76) 



As usual, the solution is parameterized by the dimensionless size of the 
plastic zone, Ap. For a given value of Ap, the corresponding load parame- 
ter is given by ( 1761) and the distribution of plastic strain is obtained by 
substituting Ci given by (1751) and C2 = into ( 167|) . 
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4-2. Solution Based on Variational Formulation 

For the variational formulation, the yield condition is used in the form 
fl2T|) instead of fl60|) . and equation fl6T]) is replaced by 



e 



-«^/^^^ (e«^/^^?«:„(0')' + ^niO = 1 - 0e-^^/^^? (77) 



As shown in Appendix D, the general solution of the corresponding homoge- 
neous equation is a linear combination of functions 

KmK) = e^,F,(^;-;l^j (78) 



wher e iFi deno t es the so-called confluent hypergeometric function of the first 



kind ISneddonI (119561 ) . defined by formulae (lD.12p - (lD.13p in Appendix D. 
A particular solution of the non-homogeneous equation (1771) could be con- 
structed by variation of constants. Following the same procedure as in Sec- 
tion 14.11 we obtain a particular solution 

UO = CiiO^nAO + C2{0^nA0 (80) 

with 

CliO = [ , J:-^^^^ . {l-<Pe-^) de (81) 

m) = - [ , , (i-0e'^) de (82) 

Unfortunately, these integrals cannot be evaluated analytically. 

An alternative approach can be based on the Green function of the dif- 
ferential operator on the left-hand side of (177j) . Formally, the Green function 
represents the solution of the differential equation with the right-hand side 
replaced by Dirac distribution centered at a point rj. At points ^ different 
from rj, the right-hand side is zero and the solution is a linear combination 
of functions (178|) - (179l) . Due to the singularity at ^ = 77, different coefficients 
of linear combination must be used for ^ < 77 and for ^ > 77. Moreover, these 
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coefficients depend on the specific value of rj. Tlierefore, we can write tlie 
Green function as 



^^^'"^^ \ i?i(r/)K„,i(0 + B^ivMO for r/ < e < ^^""^ 

and impose at ^ = 77 the continuity condition for the value and the unit jump 
condition for the first derivative. This leads to two equations, 

Ai(r/)K„,i(r/) + A2(r/)K„,2(^) = 5i(r/)/€„,i(r/) + 52(r/)K„,2(^) (84) 
A.ivKiiv) + A2ivK,2iv) = fii(r/)<i(r/) + i?2(r/X2(^)-l (85) 

In usual problems solved on a fixed interval, the Green function should also 
satisfy two boundary conditions (one at each boundary point). However, in 
our case the exact position of the boundary points is not specified and the 
number of conditions to be satisfied at each boundary point is two (vanishing 
value and vanishing first derivative). Making use of symmetry, we can restrict 
attention to non-negative values of ^ and rj and impose only one condition 
of vanishing derivative at ^ = 0, while at ^ = Ap we still need to satisfy 
two conditions. One of them can be incorporated into the Green function, 
and the other will be imposed aposteriori on the resulting solution and will 
provide a link between the load parameter and the dimensionless size of 
the plastic zone, A^. 

Based on the foregoing discussion, we constrain the Green function by 
conditions of vanishing derivative at ^ = and vanishing value at ^ = Ap, 
which gives two additional equations, 

A(^)<i(0) + ^(^)<2(0) = (86) 
Si(r/)/s:„„i(Ap) + 52(r/)K„„2(Ap) = (87) 

Since k^i(0) = (see Appendix D) and ^'^2(0) 7^ O5 equation fl86l) gives 
^2(^7) = 0. Combining equations flM]l -f l5^ and flSTjl . we can express 

. / N Kn,2{v)f^n,l{K) - ^n,l(^)^n,2(Ap) , . 

MV) - (88) 

= -^^^^^^^^ (89) 
B2iv) = ""^^'^^1^";-'^^"^ (90) 
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with the auxihary function D given by 



Div) = ['^n,l(^X,2(^) - <,l(^)«^n„2(^)] l^n,liK) (91) 

and construct the Green function by substituting this back into (1831) . The 
solution of f l77|l is then given by 

= G(e, v)dv-<P r G(e, v) e-'''/'"^ dry (92) 

JQ Jo 

This solution always satisfies conditions k^(0) = and K„(Ap) = 0, which 
have been incorporated into the Green function. The remaining condition 
to be satisfied, K'„(Ap) = 0, provides a link between the load parameter and 
the size of the plastic zone. As usual, it is more convenient to express the 
load parameter in terms of the dimensionless plastic zone size Ap than vice 
versa. Indeed, we can formally write 



<(0 = r ^'(^' r G'{^, rj) e-^'/'''^ dry (93) 

Jo Jo 

where, for simplicity, 

r'(^ „^ = ^^(^-^) = / MvXAO for < e < 

^^'^^ I 5i(^)<i(0 + fi2(^)<,2(0 for < e < Ap 

(94) 

denotes the partial derivative of the Green function G with respect to its first 
argument. From condition K^(Ap) = we obtain 

A„, A„ = -^-^ ^ ^' '[ ' 95 

4-3. Results and Discussion 

For illustration, the solution has been evaluated and plotted for a range 
of values of parameter A^. Fig. |9] indicates that the plastic zone evolves 
continuously from the weakest section and its size grows monotonically, which 
confirms admissibility of the solution. Evolution of the plastic zone profile is 
shown in Fig. [10] for A^, = 1 and 10, and the plastic part of the normalized 
load-displacement diagram is in Fig. [TTJ 
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Figure 9: Relation between load parameter and plastic zone size for a smooth (Gaussian) 
stress distribution: (a) = 1.15, 2, 10, (b) Xg = 0.95, 1.01 (from top to bottom) 



In all these figures, solid curves correspond to the variational formula- 
tion from Section 14.21 and dotted curves to the standard one from Section 
14.11 For large values of Xg, e.g. 10, both formulations give almost identical 
results. Initially, the plastic zone quickly expands and the load parameter 
increases only slightly from its value 1 at the onset of plastic yielding. The 
softening part of the load-displacement diagram is almost linear and as the 
load approaches zero, the plastic zone size tends to 27r/ (for the standard 
formulation) or to a value very close to 27r/ (for the variational formulation). 
For intermediate values of Xg, e.g. 2 or 1.15, the hardening due to structural 
effects is more pronounced for the variational formulation than for the stan- 
dard one, and the load-displacement diagram becomes more ductile. Finally, 
for small values of Xg, e.g. 1.01 or 0.95, the variational formulation gives 
unlimited hardening and unlimited expansion of the plastic zone, while the 
standard formulation still leads to limited hardening up to a finite peak load, 
followed by softening and expansion of the plastic zone to its maximum size 
27r/. 

The strong hardening effect of the variational formulation for small values 
of Xg can be attributed to the fact that the stabilizing second term in ( 12T|) 
can be expanded into HPA{x)k"{x) + HPA'{x)k'{x). The second part, de- 
pendent on the derivative of the sectional area, is neglected by the standard 
formulation but taken into account by the variational one. For the problem 
considered here, the product A'{x)k'{x) is always negative (because A{x) is 
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dimensionless coordinate dimensionless coordinate 



Figure 10: Evolution of plastic strain profile for a smooth (Gaussian) stress distribution: 
(a) X, = 1, (b) A, = 10 




normalized plastic elongation normalized plastic elongation 



Figure 11: Plastic part of normalized load-displacement diagram for a smooth (Gaussian) 
stress distribution: (a) Xg = 1.15, 2, 10, (b) Xg = 0.95, 1.01 (from top to bottom) 
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decreasing in the left half of the plastic zone and increasing in the right half, 
while K,{x) does the opposite), and so the additional term has a similar effect 
as negative k"{x). Such a term increases the yield force and thus causes 
structural hardening accompanied by an expansion of the plastic zone. 

It is important to note that cases shown in Fig.[9]D and Fig.lTTb correspond 
to an extreme nonuniformity of sectional area. Indeed, = 1 means that 
the sectional area at point x = irl, which would be the right boundary of 
the plastic zone, is e'^ ~ 139 times the area of the weakest section at the 
center of the plastic zone. So this case is only of academic interest, and is 
included here for completeness. 

5. Summary and Conclusions 

In this paper, localization of plastic strain induced by a negative plas- 
tic modulus has been studied using a variational formulation of a gradient- 
enriched plasticity model. The main points can be summarized as follows: 

1. Mathematical description of one-dimensional gradient plasticity has 
been derived using a consistent variational approach, which naturally 
provides not only the differential equation that represents the yield 
condition and needs to be satisfied inside the plastic zone, but also the 
appropriate form of the boundary and jump conditions. 

2. Taking into account the jump conditions following from a variational 
principle, a prototype problem with discontinuous data (in this spe- 
cific case, with discontinuous distribution of the sectional area) has 
been handled successfully. It has been shown that the problem has a 
physically reasonable solution, which can be constructed analytically. 

3. Two additional prototype problems, one with continuous but non- 
smooth data and the other with smooth data, have been analyzed, 
and analytical or semi-analytical solutions have been provided. The 
results have been compared to an alternative model that does not have 
a variational format. 

4. The influence of various parameters on the evolution of the plastic 
zone, on the shape of the plastic strain profile and on the resulting 
load-displacement diagram has been investigated for all three cases 
mentioned above. It has been shown that the plastic zone evolves 
from the weakest section and monotonically expands, which is in most 
cases initially accompanied by an increase of the axial force over its 
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elastic limit value. This means that the structural response exhibits 
hardening despite the softening character of the material model, which 
is related to the expansion of the yielding process into stronger parts 
of the structure induced by the gradient enrichment. 
5. The variational approach adopted here is based on the condition of non- 
negative first variation of a certain energetic functional. States that 
satisfy this condition are considered as valid solutions of the problem. 
It is expected that stable solutions correspond to local minima of the 
functional. A rigorous analysis of the second variation, leading to an 
explicit stability criterion, is presented in Appendix A for the simplest 
case of a bar with uniform properties. 

Finally, let us note that the approach elaborated here for the second-order 
gradient model and for variable cross section can be extended to the fourth- 
order gradient model and to variable material properties (such as the yield 
stress or plastic modulus). The latter extension could be useful in studies 
of the development of plastic zone near the interface of two materials with 
different properties. 

Appendix A. Second Variation and Stability Conditions 

All analyses presented in the main body of the paper have been based 
on the condition of non-negative first variation of functional 11. Of course, 
this is only a necessary but not a sufficient condition for a local minimum. 
Intuitively it can be expected that solutions corresponding to a non-negative 
first variation but not to a local minimum are unstable. Therefore, it is useful 
to investigate in more detail the changes of 11 in the immediate neighborhood 
of a "point" [u, k) that represents a valid solution, i.e., leads to a non-negative 
first variation 611 for all admissible variations 6u and 6k. 

If 611 > for given 6u and 6k,, the value of 11 increases at least locally (for 
sufficiently small variations), and the second variation does not need to be 
evaluated. However, for those 6u and 6k, that give a vanishing first variation 
611, the sign of the increment is determined by the second variation. So, 
as a first step, we restrict attention to those combinations 6u and 6k for 
which 611 = 0. Since the solution {u, k) satisfies the equilibrium condition 
f lTSj) . which is an equality, the first integral in f|T6|) vanishes for an arbitrary 
displacement variation 6u. On the other hand, the yield condition f fT9|) is 
an equality satisfied in the plastic zone Ip only. Outside the plastic zone. 
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the plastic admissibility condition (12U]) becomes an inequality. Leaving aside 
some degenerated cases of neutral loading, this condition is typically satified 
as a strict inequality. To make sure that the first variation 611 vanishes, the 
plastic strain variation 6k must be set to zero outside the plastic zone. By 
similar arguments based on the jump terms, it can be shown that the values 
of 6k, and of its derivative on the boundary of the plastic zone must vanish. 

Since functional 11 given by f[T^ is quadratic, its second variation is easily 
expressed as 

6^U{6u,6k;u,k) = j \EA{6u' -6Kf dx + j \HA{6k^-1^6k'^) dx (A.l) 

Taking into account that 6k vanishes outside Xp, as justified above, the inte- 
gration domains can be reduced. After simple rearrangements (with moduli 
E and H considered as constants), the resulting stability condition can be 
written as 

I A6u'^dx+ I A{6u' -6Kf dx-'^ I A{1^6k''' -6k^) da; > (A.2) 

In a loading test performed under displacement control, the displacements 
on the physical boundary dL are prescribed and the variations 6u vanish on 
dC Condition (]A.2p should be satisfied for all such 6u and for all variations 
6k with vanishing values and vanishing first derivatives on the elasto-plastic 
boundary dXp. 

General analysis of condition (1A.2I1 for variable sectional area A would 
be very tedious. However, the case of a uniform bar (with A = const.) is 
manageable, because flA.2l) simplifies to 

f 6u'^ dx+ f {6u' - 6k)^ [ (^''^'^" - '^'^') ^ (^-3) 

In this case, the plastic zone Xp is an interval of length Lp = 2tiI (see Sec- 
tion [2]T]), and the elastically unloading zone £\Xp is a union of two intervals 
of total length L — Lp. Individual integrals in (1A.3P can be estimated from be- 
low. For the second integral, we can exploit the Cauchy-Schwarz inequality, 
which implies that 

f {6u' - 6Kf dx > ^ I / {6u' - 6k) dx ) = Lp{6e - 6k)^ (A.4) 
Jxp Lp \Jxp J 
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where 



6u' dx, 



5k, = — 



Sk dx 



(A.5) 



are constants that represent the mean values of 6u' and 6k, over the plastic 
zone. The first integral in f lA.Sp can be estimated in a similar fashion, taking 
into account that the mean value of 6u' in C\Xp is directly related to 6e, 
because of the compatibility constraint 6u' dx = 0: 



6u dx 



6u dx 



6u'^ dx > 



L-Lr 



-Lp5e 



Su' dx 



^2 



L — Lr, 



(A.6) 



(A.7) 



Finally, the last integral in flA.3|) can be estimated using the Wirtinger in- 
equality, 

which is a one-dimensional case of the Poincare inequality, with explicitly 
known optimal constant. The theorem is applicable to the zero- mean part 
of 5k and implies that 

I 6k'^ dx> I {5k - 6kY dx= I 6k^ dx - Lp5R^ (A.8) 

Since Lpjlix = I, we get 

[ {P6k'^ - 6k'^) dx > -Lp6K^ (A.9) 



Substituting flA.4p . flA.7p and (IA.9P into (lA.Sp . we obtain condition 



H 



L-L, 



6e' + Lp{6e - dRy + —LpdR' > 

E 



(A.IO) 



which means that a certain quadratic form of variables Se and Sk should be 
positive semidefinite. This leads to the following restrictions on the param- 
eters: 



r2 



L-L, 



+ Lp>0 
H 



L> Lr 



(A.ll) 



Lp + —Lp>0 E + H>0 (A.12) 
E 



^P 



L-Lr 



LP + ^-^p 



^>-f (A,13) 
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The first restriction corresponds to our tacit assumption that bar is longer 
than Lp, so that the full plastic zone can develop (for shorter bars, the anal- 
ysis would have to be modified). The second restriction excludes snapback 
of the stress-strain diagram with no gradient effects. The third restriction 
is the most stringent one. It guarantees stability of the localized solution 
under displacement control and can be interpreted e.g. as a constraint on 
the maximum length of the bar (with respect to the characteristic length /, 
reflected by the plastic zone size Lp = 27r/). It is reassuring that this condi- 
tion exactly coincides with the condition of negative slope of the post-peak 
load-displacement diagram. Indeed, the total elongation of the bar after the 
onset of plastic yielding can be expressed as a sum of the elastic and plastic 
parts: 

f f F f FL 

Utot= I s dx = I — — dx + / K dx = —— + Up (A.14) 



c J c 



EA r EA 



Substituting the plastic strain distribution according to which refers to 
the plastic zone Xp = (— ttZ, nl), the plastic part of elongation turns out to be 

. r (l ^ CO. f ) dx ^ ^2.1 ^ ,A.15) 

J-TVI 



H \ U H HA H 



The slope of the post-peak part of load-displacement diagram is the reciprocal 
value of the tangent structural compliance L/EA + Lp/HA. Snapback occurs 
if the tangent compliance (and thus also the tangent stiffness) is positive, i.e., 
if 

| + |>0 (A.16) 

This is of course equivalent to L / Lp < —H/E, which holds if and only if the 
stability condition (IA.13|1 is violated. 

To illustrate the difference between the condition of non-negative first 
variation and the (stronger) condition of a local minimum, let us recall that 
the localization problem for a uniform bar admits not only solutions with 
the plastic zone of size Lp = 27il, but also solutions with Lp equal to integer 
multiples of 27r/; see Section 12.11 For instance, a plastic strain distribution 
given by 

K{x) = i ^i^(^-"°'T) forxGX, = (-27r/,2vr/) ^^^^^^ 
for a; ^ Xp 
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satisfies differential equation fl23p as well as conditions k = and k' = 
at the boundary of Xp, i.e., at points ±27r/. The length of the plastic zone, 
Lp = Airl, is now the double of the minimum plastic zone length. Consider 
an admissible variation of k in the form 

6k{x) = I ' ^Sn(x) (l - cos ^) for x G 1^ ^^^^^^ 
[ for X ^ Xp 

where c is an arbitrary constant, not exceeding in magnitude the positive 
constant {a — Cq)/!! (for larger |c|, the function n + 6k would not be non- 
negative and would not belong to the space of admissible functions K; defined 
in (fT4|) ). A simple calculation reveals that 

/ {1^6k'^ - 6k^) dx = 

Jlp J ~2-kI 

(A.19) 

Since the mean value of 5k is zero, the displacement variation 5u can be 
selected such that 5u' = 5k, and then the first two integrals in flA.31) vanish. 
The third integral, evaluated in flA.191) . is negative and is multiplied by a 
positive constant —H/E, and so the stability condition ( ]A.3p is always vio- 
lated. This proves that solutions with plastic zone exceeding the minimum 
length 2ttI would be unstable, independently of the total bar length L. 



X 



sm 



X 



cos ■ 



dx = —4c ttI 



Appendix B. General Solution of Equation ( 1461) 

For negative ^, equation (j4]) can be written as 

Using substitutions 

^ = V-K (B-2) 

k{0 = {Xg + 09i\ + = r]9{v) (B.3) 

where g is a new unknown function and 77 is a shifted dimensionless spatial 
coordinate, we can convert (IB.ip into 

V^9"iv) + Vg'iv) + - ^)9iv) = V- (B.4) 

^9 
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where primes denote derivatives with respect to rj. The homogeneous coun- 
terpart of (IB.4p is the Bessel equation of order 1^ = 1, and its general solution 
can be written as 

9hiv) = CMr]) + C2Y^{r]) (B.5) 

where Ci and C2 are arbitrary constants, and Ji and Yi a r e resp ectively the 
Bessel functions of the first and second kind; see iKorenevI (120021 ) . 

Now we need to find a particular solution for the given right-hand side. It 
turns out that the expression on the left-hand side of (IB .41) gives a multiple 
of ?7 if is set simp ly to l/?7 , and a multiple of 77^ if g is set to the Struve 
function Hi(x); see 
solution in the form 



KorenevI ( 120021 ) . Therefore, we look for the particular 



ki 

giji) = h k2Ui{v) 

V 



(B.6) 



and substituting into the left-hand side of (]B.4|) we get the condition 



kiT] + 



2k2r]^ 



rj- —T] 



(B.7) 



from which ki = 1 and k2 = — 7r0/2Ag. Combining the particular solution g 
given by ( IB. 61) with the general solution of the homogeneous equation gh given 
by ( IB. 51) and substituting this back into ( IB.Sp . we get the general solution of 
( ]B.1|) in the form 

<0 = (A,+0[C^iJi(A,+0+C2Yi(A,+0] + l-^(A,+0Hi(A,+0 (B.8) 



Recall that (IB.ip is the specific form of (|4]) valid for ^ < 0. For ^ > 0, it is 
sufficient to replace A^ + ^ by Aj, — ^. The resulting expression valid for both 
positive and negative ^ is given in (l47l) . 



Appendix C. General Solution of Equation ( 1611) 

The integrals in (IBS]) - (|M]) are most conveniently evaluated if one intro- 
duces an auxiliary complex function 



c{0 = Ci{0-^C2{0 



(C.l) 
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Substituting from one gets 

= cos^ — 2 sin^ + J e 



^'/2^9(sin^ + 2cosO 



e '^ + ic 



Since 



2 



the last integral in (]C.2p can be written as 



exp 



(C.2) 



(C.3) 



de (C.4) 



and substitution ^ = \/2Xgt — iXl leads to 



exp 



AoW— erf, ^ 



3 



(C.5) 



where erf is the "error function" defined by 



erf(x) 



(C,6) 



Combining ( ]C.4[) and ( IC.Sj) and substituting back into (1C.2I) . we obtain 



(7(O = e-^ + .0A,^e-^/^erf(^ 



(C.7) 



Functions Ci and C2 could now be extracted from the real and imaginary 
part of the expression in flC.7|) . But this is not even necessary, because the 
particular solution of equation f l6T|) given by formula (|62|) can be presented 

as 



k^iO = Re (^1(0 - zC'2(0)(cose + zsinO 



Re 



C(e)e^« (C.8) 
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where Re stands for the real part. Substituting for C{^) according to (lC.7p . 
we get the particular solution in the form 



^n(0 = 1 + 



-"^/2Re 



,erf( ^i^]e^« 



V2\ 



Recalling the definition of the Dawson function lOlverl (jl997( ) 



(C.9) 



Fix) = e-"' / e''dt = -!^e-"'erf(zx) 



(C.IO) 



and taking into account that F{x) = F{x) (with the bar denoting the com- 
plex conjugate), we can rewrite the result as 



+ F 



(C.ll) 



The general solution ( 167|) of equation (!6T]) is then obtained by adding a linear 
combination of functions cos^ and sin^ with arbitrary coefficients. 



Appendix D. General Solution of Homogeneous Form of Equa- 
tion ( 177]) 



The homogeneous counterpart of equation (1771) can be written as 



(D.l) 



Defining a rescaled spatial variable r] = ^/Xg\/2 and expressing the unknown 
function as 



«n(0 = e ''''ou 



XaV2 



e u{rf) 



(D.2) 



where u is a transformed unknown function, we can convert fID.ll) into the 
so-called Hermite differential equation. 



m"(?7) - 2?7m'(?7) + 2vu{'q) = 



(D.3) 



where z/ = — 1 is a real parameter (in our case larger than —1) and primes 
denote differentiation with respect to t]. 
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The solution of the Hermite equation (1D.3P can be constructed in terms 
of infinite power series 



r=0 

Expressing the derivatives 

oo 

r=l 

oo 

u"{7]) = ^a,s(s- l)r/^-2 (D.6) 

s=2 

and substituting them into (ID.Sp . we obtain 

oo oo oo 

ass{s - l)r]'-^ - 2 XI + ^ arv"' = (D.7) 

s=2 r=l r=0 

In the first sum, s can be replaced by r + 2, with r running from to infinity, 
and in the second sum, r can also run from without changing the result 
(because of the presence of the factor r which anihilates the term with r = 0). 
The equation can thus be rewritten as 

oo 

[ar+2{r + 2)(r + 1) - 2a,,r + 2z/a,.] r]'' = (D.8) 

r=0 

and the term in the square brackets must vanish for each individual value of 
r = 0,l,2,.... This condition results into the recursive formula 

2(r-u) 
(r + l)(r + 2) 

Note that ar+2 is expressed in terms of a^. It is therefore possible to select 
arbitrary values of Oq and ai, and then express all other coefficients with 
even subscripts in terms of qq and all other coefficients with odd subscripts 
in terms of ai. Setting oq = 1 and ai = 0, or oq = and ai = 1, leads to two 
linearly independent functions 

2z/ 5 2V(z/-2) , 2V(z/-2)(i/-4) . 

(D.IO) 

2(z/-l) , 22(iy-l)(z/-3) , 

Mv) = ^-^^^ +^ — +••■ 
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and every solution of equation (1D.3P can be expressed as their linear combi- 
nation. The fundamental solutions can be conveniently expressed in terms 
of the so-called confluent hypergeometric function of the first kind, denoted 
as 7; x) (note that the symbol F has a left subs c ript a nd a right sub- 

script), which is defined by the infinite series ISneddonI (Il956l ) 



iFi(a;7;x) = ^ 



r=0 



(D.12) 



Here, {•)r is the so-called Pochhammer symbol, which can be expressed in 
terms of Euler's gamma function: 



(«)r 



a(a + 11 



(a 



r{a + r) 

r(a) 



(D.13) 



It is easy to verify by simple substitution that the fundamental solutions ui 
and U2 from ( 1D.10l) -( lD.lll) can be rewritten as 



uiiv) = iFi(-z//2;l/2;r/2) 
U2iv) = r/iFi((l-z/)/2;3/2;r/2) 



(D.14) 
(D.15) 



Substituting this into (1D.2P and replacing z/ by — 1 and 77 by ^/\g\/2, we 
finally obtain two linearly independent solutions of f ID.ip in the form 



«n,l(0 
«n,2(0 



iFl 

V2\ 



-iFi 



: '2'2A2 



2' 2A2 



(D.16) 
(D.17) 



Note that function ui is even and U2 is odd, and so K,n,i is even and Kn,2 is 
odd. One useful consequence is that i(0) = 0. 
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